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Abstract: A hybrid experiment has been started by AS7 collaboration at Tibet, China, since May 2009, that consists of 
a burst-detector-grid (YAC, Yangbajing Air shower Core array) and the Tibet air-shower array (Tibet-Ill). The first step 
of YAC, called YAC-I, contains 16 detector units and observes high energy electromagnetic particles in air-shower cores 
within several meters from the shower axis, and Tibet-Ill array measures the total energy and the arrival direction of air 
showers. YAC-I is used to check hadronic interaction models currently used for air-shower simulations such as QGSJET, 
SIBYLL , EPOS etc. through the multi-parameter measurement in air-shower cores. In this paper, we used a data set 
collected from May 1st 2009 through February 23rd 2010 by the YAC-I. The effective live time used for the present 
analysis is 169.65 days. The preliminary results of the interaction model checking at *10 TeV energy region is reported 
in this paper. 



Keywords: air shower core, hadronic interaction model,cosmic ray 



J.SHAO et al. TEST OF THE HADRONIC INTERACTION MODELS AROUND * 10 TVV 



1 Introduction 

The interpretation of the extensive air showers (EAS) is 
known to inevitably depend on the Monte Carlo simula- 
tions which are based on some hadronic interaction mod- 
els and cosmic ray composition models. At present, the 
simulation code CORSIKA that is comprehensively used 
in the surface cosmic ray studies includes many interaction 
models. For multi-parameter measurements of EASs, it is 
known that no interaction model can explain all data con- 
sistently. Therefore, the hadronic interaction models need 
to be checked and further improved. 

In this paper, we report our approach and the preliminary 
results to check the hadronic interaction models QGSJET2 
and SIBYLL2.1 at an energy region of *10 TeV using the 
data obtained by the newly constructed YAC-I (Yangbajing 
Air shower Core detector, the first stage). The energy re- 
gion of *10 TeV is chosen by the following considerations: 
1) The primary composition at this energies has been better 
measured by direct measurements, and the uncertainty[l] 
is smaller. 2) The corresponding energies in the center-of- 
mass system are around *100 GeV, for that SppS collider 
made good measurements on the inelastic cross section and 
on the particle production in the near-forward region. The 
uncertainty from the extrapolation to the very forward re- 
gion is relatively small. 3) The check of interaction mod- 
els should be step-by-step executed from lower energies to 
higher energies. 4) Due to the high altitude of our Tibet 
experiment we could have good EAS core events recorded 
by YAC-I at *10 TeV region. 

2 YAC-I experiment 

YAC-I, consisting of 16 scintillation detectors, has been 
successfully operated at Yangbajing in Tibet, China, since 
May 2009, together with the Tibet-Ill EAS array. Each unit 
of YAC-I is composed by a lead layer of 3.5 cm thickness 
(~7 r.l.) and a plastic scintillator of size 40 cm x50 cmx 1 
cm. Each YAC-I detector is used to record the electro- 
magnetic showers induced by high energy electrons and/or 
photons in the EAS cores. 16 detectors are arranged by 
4x4, forming a ~ 10 m 2 covering. To achieve a wider 
burst size measurement under the lead layer, a wide dy- 
namic range from 1 MIP (Minimum Ionization Particle) to 
10 6 MIPs is demanded which is realized by two photomul- 
tipliers(PMT), i.e. a high-gain PMT and a low-gain PMT 
for the range of 1 ~ 3 x 10 3 and 10 3 ~ 10 6 MIPs, respec- 
tively. The response linearity of each YAC-I detector was 
calibrated by cosmic-ray single muons and by the accel- 
erator beam (BEPC-LINAC)[10]. YAC-I is triggered when 
any one of 16 detectors records a local shower with the size 
of at least 20 MIPs. The event rate is about 30 Hz. The total 
live time of our data set in present analysis is 169.65 days. 

3 Simulations 

A Monte Carlo simulation has been carried out on the de- 
velopment of EASs in the atmosphere and their response in 
YAC-I. The simulation code CORSIKA (version 6.204)[2] 
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Figure 1: Primary cosmic -ray composition for the HD 
model and the NLA model. The all-particle spectrum, 
which is the the sum of all components, is nomalized to 
the Tibet data and compared with other experiments. Black 
solid lines denoted p and He spectrum that are fitted from 
PAMELA and CREAM data [6] [8] [9]. 

including QGS JET2 and SIBYLL2. 1 hadronic interaction 
models are used to generate air shower events. In this 
work, two primary composition models are used, the heavy 
dominant model(HD)[3] and the non-linear acceleration 
model(NLA)[4], as shown in Fig.l. The proton spectrum 
of the two models is connected with the direct experiment 
in the low energy and consistent with the spectrum obtained 
from the Tibet AS+EC experiment in the high energy. The 
He spectrum of HD model coincides with the results from 
RUNJOB and ATIC-I, but the He spectrum of NLA co- 
incides with the results from JACEE, ATIC-II, CREAM3. 
The sum of all single-component spectra can reproduce the 
sharp knee in all particle spectrum [4]. The fractions of the 
component of the two composition models in different en- 
ergy regions are listed in Table 1 . 

In the simulation, the primary energy is sampled from 1 
TeV to infinite with zenith angles from to 60 degrees inci- 
dent isotropically. The axis of each EAS event is randomly 
dropped onto an area of 32.84 m x32.14 m with YAC-I at 
its central part. When high energy electrons or photons hit 
a YAC-I detector, the EPICS(8.64)[5] code is used to gener- 
ate the cascade showers in the Pb layer and in the detector. 
To identify an AS core event we use following quantities: 

Nf,: number of shower particles recorded by a YAC-I de- 
tector; 

Nhit- number of YAC-I detectors with Nf, higher than a 
threshold value Nf, m i n (If Nf,>Nf, m in for a YAC detector, 
it is called that this detector is fired); 

"^Nf,: the sum of all Nf, from 1 to Nhit', 

Nl° p : the maximum of all Nf,; 

< R >: the mean lateral distance from the center of a fired 
detector to the Nf, weighter center of all fired detectors; 

< Nf,xR >: the mean Nf, weighted lateral distance from 
the center of a fired detector to the Nf, weighted center. 

By choosing a Nf, m i n value and setting some event selec- 
tion conditions we can obtain different event samples for 
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Composition model 


Component 


1 - 10 TeV 


10 - 100 TeV 


10 3 -10 4 TeV 




P 


38.6% 


32.0% 


24.2% 


HD 


He 


24.7% 


22.4% 


19.1% 


model 


M 


24.7% 


27.1% 


27.3% 




Fe 


10.4% 


18.5% 


29.4% 




P 


47.3% 


31.1% 


26.3% 


NLA 


He 


30.2% 


25.1% 


28.9% 


model 


M 


30.3% 


32.4% 


34.4% 




Fe 


12.8% 


11.3% 


10.6% 



Table 1 : The fractions of the components of the HD model and the NLA model. 



Mode energy 


- 35 TeV 


- 70 TeV 


- 90 TeV 


QGSJET2+HD 


2640 


8893 


5352 


SIBYLL2.1+HD 


2921 


9239 


5495 


QGSJET2+NLA 


4484 


15515 


9312 


SIB YLL2 . 1 +NL A 


3941 


12568 


7492 


YAC-I data 


1773 


5640 


3220 



Table 2: The numbers of core events of three samples selected by conditions 1), 2), 3) (see Section 3)in Monte Carlo 
simulation and YAC-I data. 



that primary energies range at different region. For the 
present work using Nb m i n =20Q and the conditions of 1) 
N hit =3, J2N b > 3500, 2) N hit =4, > 1700, 3) 

N hit =5, Y, N b > 2200, three Monte Carlo samples with 
the mode energy at -35 TeV, -70 TeV and -90 TeV are 
obtained, respectively. Their sample sizes are seen from 
Table 2. In all 12 cases, Monte Carlo shows that the core 
resolution is better than 2 m if the Nf, weighted center is 
taking as the AS core. 

4 Data analysis 

The present analysis uses a subset of data collected in the 
period from May 1st, 2009 to February 23rd, 2010, having 
stable instrument conditions. We first check some noises 
appeared in our data set due to the hardware and the en- 
vironment conditions. It is found that most noises appear 
as smaller 'signals'. Taking A^, mirl =200 can remove them 
and achieve a biggest available data sample. In addition, we 
found some gain drift for some YAC-I detectors during the 
operating process, by checking Nb spectrum of each YAC- 
I detector. By an 'off-line self-calibration' method this ef- 
fect is carefully treated and corrected. Then the fit between 
high-gain signals and low-gain signals in their overlapping 
region is executed and N h (or MIPs) is obtained for each 
fired YAC-I detector. 

After the off-line calibration, — 150000 events with Nb > 
200 and Nhit > 1 are obtained. Three experimental data 
samples with the same selection criteria of Monte Carlo 
are obtained. The sizes of each sample are also listed in 
Table2. 



5 Results and Discussion 

Since our Monte Carlo simulation is started from 1 TeV, 
in order to normalize MC data and experimental data, we 
need to know the integral intensity of all particles of cos- 
mic rays at E = 1 TeV. Starting from Horandal's spectra of 
each composition[7], we improve the major 8 ones (p, He, 
C, O, Ne, Mg, Si , and Fe) by the newest measurements 
[8][9][10]. The resultant integral intensity: I(E >1 TeV) = 
0.139 cm-W- 1 with the error +0.0013, -0.0012 com- 
ing from the error of the index of each of the 8 spectra. 

The comparison of our data with Monte Carlo simulations 
of QGSJET2+NLA, QGSJET2+HD, SIBYLL2.1+NLA 
and SIBYLL2.1+HD is seen from Fig.2 and Fig.3: 

(1) The shape of the distributions of Nbl^Nb are consis- 
tent between the YAC-I data and simulation data in all four 
cases, indicating that in the * 10 TeV energy region the par- 
ticle production spectrum of QGSJET2 and SIBYLL2. 1 
may correctly reflect the reality within our experimental 
systematic uncertainty of a level about 10%. 

(2) The comparison of event absolute intensities in all 
cases, as seen from Fig.2 and Fig.3, shows some discrep- 
ancies. The smallest one is SIBYLL2.1+NLA and the 
most obvious one is QGSJET2+HD. For a given interaction 
model the NLA composition is better than HD. For a given 
composition model SIBYLL2.1 is better than QGSJET2. 

(3) But note that, as seen from Fig 1, both composition 
models NLA and HD used a steeper He spectrum, if com- 
paring with the new results from PAMELA and CREAM. 
It can be estimated that NLA under-estimated the number 
of He events for about 30-40% in the related energy region. 
If involving this factor, QGSJET2 results can go higher and 
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Figure 2: The absolute intensity of Nbf^Nb distribu- 
tion of the three samples for QGSJET2+HD model, 
SIBYLL2.1+HD model, QGSJET2+NLA model, 
SIBYLL2.1+NLA model and the experimental data. 

be possible better consistent with data. A further analysis 
is going on. 

6 Summary 

A smaller high-altitude AS core detector YAC-I shows the 
ability and sensitivity in checking the hadronic interaction 
models in *10 TeV region. 

The experimental distributions, AVX^h an ^ nas tne 
shape very close to the Monte Carlo predictions of 
QGSJET2+NLA, QGSJET2+HD, SIBYLL2.1+NLA and 
SIBYLL2.1+HD. Some other quantities, such as N b , J2 N b, 
Nb tov , Rw, RwxJVt have the same behavior as well, 
though we did not show them in this paper due to the limit 
of the space. 

Some discrepancies in the absolute intensities are seen. 
Data normally shows a higher intensity than Monte Carlo. 
Taking a more hard He spectrum (and somewhat hard pro- 
ton spectrum) as given by CREAM at the 1 - 1 00 TeV region 
can improve this situation. A further study is going on. 



Figure 3: The ratio of absolute intensity between MC and 
Exp. data. 
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